library(data.table)
library(magrittr)
library(ggplot2)
theme_set(theme_minimal())

metadata <- fread("ihmp/data/ibdmdb_healthy.csv")
sra_data <- fread("ihmp/data/ibdmd_healthy_runtable.tsv")[, .(Run, `Library Name`)]
sra_data[, "External ID" := tstrsplit(`Library Name`, "_MGX")[[1]]]
metadata <- metadata[sra_data, on="External ID", nomatch=0]
metadata
man <- fread("db/data/manifest.csv")[, .(db, rank, matched_taxid, seqlength)]

foods <- fread("ihmp/data/food_abundance.csv")
foods <- man[foods, on="matched_taxid"]
foods[, "tpm" := 1.0e6 * reads / as.double(seqlength)]
foods[is.na(tpm), tpm := 0.0]

contents <- fread("ihmp/data/food_content.csv")
contents
tests <- list(
  vegetables = list(group="Vegetables (salad, tomatoes, onions, greens, carrots, peppers, green beans, etc)", ex=expression(food_group == "Vegetables")),
  fruits = list(group="Fruits (no juice) (Apples, raisins, bananas, oranges, strawberries, blueberries", ex=expression(food_group == "Fruits")),
  beans = list(group="Beans (tofu, soy, soy burgers, lentils, Mexican beans, lima beans etc)", ex=expression(food_subgroup == "Beans")),
  `white meat` = list(group="White meat (chicken, turkey, etc.)", ex=expression(food_subgroup == "Poultry")),
  shellfish = list(group="Shellfish (shrimp, lobster, scallops, etc.)", ex=expression(food_subgroup %chin% c("Mollusks", "Crustaceans"))),
  fish = list(group="Fish (fish nuggets, breaded fish, fish cakes, salmon, tuna, etc.)", ex=expression(food_subgroup == "Fishes")),
  `red meat` = list(group="Red meat (beef, hamburger, pork, lamb)", ex=expression(food_subgroup %chin% c("Swine", "Bovines", "Caprae")))
)
library(patchwork)

tables <- list()

freqs <- c(
  `No, I did not consume these products in the last 7 days` = 0,
  `Within the past 4 to 7 days` = 1/5.5,
  `Within the past 2 to 3 days` = 1/2.5,
  `Yesterday, 1 to 2 times` = 1.5,                                
  `Yesterday, 3 or more times` = 3,
  `NA` = 0
)

results <- list()
plots <- list()
merged <- metadata[foods, on=c(`Run` = "sample_id")]
for (i in 1:length(tests)) {
  name <- names(tests)[i]
  vals <- tests[[i]]
  dt <- merged[, .(abundance=sum(reads[eval(vals$ex)]), group=.SD[[vals$group]][1], total_reads=total_reads[1]), by="External ID"]
  ns <- dt[, .N, by="group"] |> setkey(group)
  dt <- dt[group %chin% names(freqs) & ns[group, N] >= 10]
  dt[, "frequency" := freqs[group]]
  dt[, group := factor(group, levels=names(freqs)[names(freqs) %in% group])]
  dt[, "relative" := (abundance + 1) / (total_reads)]
  dt[, "item" := name]
  plots[[name]] <- ggplot(dt) + aes(x=group, y=log10(relative)) + 
    geom_jitter(width=0.3) + 
    geom_boxplot(width=0.1) + 
    stat_smooth(aes(group=1), method="lm") + 
    labs(title=name, x="consumption frequency [servings/day]", y="abundance")
  print(name)
  abmod <- lm(log10(relative) ~ frequency, data=dt)
  premod <- glm((abundance > 0) ~ frequency, data=dt, family=binomial(link="logit"))
  tables[[name]] <- dt
  results[[name]] <- dt[, .(n=.N, relative=mean(log10(relative*1e6)), sd=sd(log10(relative*1e6)), detected=sum(abundance > 0), item=name), by="group"]
  results[[name]][, "p_abundance" := summary(abmod)$coefficients[2, 4]]
  results[[name]][, "p_prevalence" := summary(premod)$coefficients[2, 4]]
}
[1] "vegetables"
[1] "fruits"
[1] "beans"
[1] "white meat"
[1] "shellfish"
[1] "fish"
[1] "red meat"
r <- rbindlist(results)
ta <- rbindlist(tables)
ggplot(r) +
  aes(x=group, y=relative) +
  geom_jitter(data=ta, aes(y=log10(relative*1e6)), stroke=0, size=0.5, width=0.3) +
  geom_pointrange(aes(ymin=relative-sd, ymax=relative+sd), shape=21, fill="white") +
  geom_text(data=unique(r, by="item"), aes(x=0, y=Inf, label=sprintf("p[lm]=%.2g\np[logit]=%.2g", p_abundance, p_prevalence)), size=3, vjust=1, hjust=0, nudge_y=0.5) +
  facet_wrap(~ item, ncol=2) +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1)) +
  labs(x="", y=expression(paste("log relative abundance [",log[10],RPM, "]"))) +
  guides(color=F)

ggsave("figures/ihmp_foods.pdf", width=2.5, height=9)
ggplot(foods[seqlength>1e6], aes(x=seqlength/1e9, y=relative_abundance)) + geom_point() + stat_smooth(method="lm", color="tomato") +
  labs(x="haploid genome size [Gbps]", y="relative food abundance") + scale_x_log10() + scale_y_log10()

foods[seqlength>1e6, cor.test(reads, seqlength)]

    Pearson's product-moment correlation

data:  reads and seqlength
t = 1.9748, df = 2209, p-value = 0.04841
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
 0.000294607 0.083521188
sample estimates:
       cor 
0.04198072 
ggsave("figures/genome_size.pdf", width=5, height=3)

Diet distances

species <- fread("ihmp/data/S_counts.csv")[grepl("Bacteria|Archaea", lineage)]
microbiome_mat <- dcast(species, sample_id ~ lineage, value.var = "new_est_reads", fill=0, fun.aggregate = sum)
sids <- microbiome_mat[, tstrsplit(sample_id, "S_")[[2]]]
microbiome_mat <- as.matrix(microbiome_mat[, "sample_id" := NULL])
rownames(microbiome_mat) <- sids
good <- colMeans(microbiome_mat) > 10
microbiome_mat <- microbiome_mat[, good]
# micro_rare <- phyloseq(otu_table(microbiome_mat, taxa_are_rows = F)) |> rarefy_even_depth(100000) |> otu_table()
microbiome_relative <- microbiome_mat / rowSums(microbiome_mat)
diet <- metadata[, names(metadata)[72:92][-11], with=F]
diet <- apply(diet, 2, function(x) ifelse(x == "", 0, freqs[x]))
rownames(diet) <- metadata[["Run"]]
diet <- diet[rowSums(diet) > 0, ]
diet <- diet[, colMeans(diet >0) > 0.1]
dim(diet)
[1] 362  19
foodab <- foods
foodab[, "id" := paste(matched_taxid, species, wikipedia_id, sep="|")]
food_mat <- dcast(foodab, sample_id ~ id, value.var = "reads", fill=0, fun.aggregate = sum)
sids <- food_mat[, sample_id]
food_mat <- as.matrix(food_mat[, sample_id := NULL])
rownames(food_mat) <- sids
food_relative <- food_mat[rowSums(food_mat) > 0, ]
food_relative <- food_relative / rowSums(food_relative)
food_relative <- food_relative[rowSums(food_relative) > 0, ]
food_relative <- food_relative[, colMeans(food_relative >0) > 0.1]
dim(food_relative)
[1] 359  16

Beta diversity tests

Let’s write a little function that runs the test for microbiome <-> other comparison and plots and returns the results.

library(vegan)

micro_mantel <- function(first, other, firstname, othername) {
  sids <- intersect(rownames(first), rownames(other))
  first_dist <- vegdist(first[sids, ], "bray")
  other_dist <- vegdist(other[sids, ], "bray")
  
  test <- mantel(first_dist, other_dist, method="pearson")
  results <- data.table(first=as.numeric(first_dist), other=as.numeric(other_dist))
  pl <- ggplot(results) +
    aes(x=first, y=other) +
    geom_point(size=1, alpha=0.25, stroke=0) +
    stat_smooth(method="lm") +
    labs(x=paste(firstname, "distance [Bray]"), y=paste(othername, "distance [Bray]"))
  print(pl)
  return(data.table(x=firstname, y=othername, r=test$statistic, p=test$signif, perms=tests$permutations))
}

Now we run it for the measures.

mantel_tests <- list(
  micro_mantel(microbiome_relative, diet, "microbiome", "FFQ"),
  micro_mantel(microbiome_relative, food_relative, "microbiome", "MEDI"),
  micro_mantel(diet, food_relative, "FFQ", "MEDI")
) |> rbindlist()

mantel_tests

Figure data

fwrite(ta, "figure_data/iHMP_points.csv")
fwrite(r, "figure_data/iHMP_summaries.csv")
fwrite(mantel_tests, "figure_data/iHMP_mantel.csv")

fwrite(foods[seqlength>1e6], "figure_data/ihmp_genome_size.csv")

Analysis for hibiscus reads

library(patchwork)

tables <- list()

freqs <- c(
  `No, I did not consume these products in the last 7 days` = 0,
  `Within the past 4 to 7 days` = 1/5.5,
  `Within the past 2 to 3 days` = 1/2.5,
  `Yesterday, 1 to 2 times` = 1.5,                                
  `Yesterday, 3 or more times` = 3,
  `NA` = 0
)
tests <- names(metadata)[72:92][-11]

results <- list()
plots <- list()
merged <- metadata[foods, on=c(`Run` = "sample_id")]
for (i in 1:length(tests)) {
  name <- tests[i]
  dt <- merged[, .(abundance=sum(reads[genus=="Hibiscus"]), group=.SD[[name]][1], total_reads=total_reads[1]), by="External ID"]
  ns <- dt[, .N, by="group"] |> setkey(group)
  dt <- dt[group %chin% names(freqs) & ns[group, N] >= 10]
  dt[, "frequency" := freqs[group]]
  dt[, group := factor(group, levels=names(freqs)[names(freqs) %in% group])]
  dt[, "relative" := (abundance + 1) / (total_reads)]
  dt[, "item" := name]
  if (dt[, uniqueN(group) == 1])  {
    next
  }
  plots[[name]] <- ggplot(dt) + aes(x=group, y=log10(relative)) + 
    geom_jitter(width=0.3) + 
    geom_boxplot(width=0.1) + 
    stat_smooth(aes(group=1), method="lm") + 
    labs(title=name, x="consumption frequency [servings/day]", y="abundance")
  print(name)
  abmod <- lm(log10(relative) ~ frequency, data=dt)
  premod <- glm((abundance > 0) ~ frequency, data=dt, family=binomial(link="logit"))
  tables[[name]] <- dt
  results[[name]] <- dt[, .(n=.N, relative=mean(log10(relative*1e6)), sd=sd(log10(relative*1e6)), detected=sum(abundance > 0), item=name), by="group"]
  results[[name]][, "p_abundance" := summary(abmod)$coefficients[2, 4]]
  results[[name]][, "p_prevalence" := summary(premod)$coefficients[2, 4]]
}
[1] "Soft drinks, tea or coffee with sugar (corn syrup, maple syrup, cane sugar, etc)"
[1] "Diet soft drinks, tea or coffee with sugar (Stevia, Equal, Splenda etc)"
[1] "Fruit juice (orange, apple, cranberry, prune etc.)"
[1] "Water"
[1] "Alcohol (beer, brandy, spirits, hard liquor, wine, aperitif, etc.)"
[1] "Yogurt or other foods containing active bacterial cultures (kefir, sauerkraut)"
[1] "Dairy (milk, cream, ice cream, cheese, cream cheese)"
[1] "Fruits (no juice) (Apples, raisins, bananas, oranges, strawberries, blueberries"
[1] "Vegetables (salad, tomatoes, onions, greens, carrots, peppers, green beans, etc)"
[1] "Beans (tofu, soy, soy burgers, lentils, Mexican beans, lima beans etc)"
[1] "Whole grains (wheat, oats, brown rice, rye, quinoa, wheat bread, wheat pasta)"
[1] "Starch (white rice, bread, pizza, potatoes, yams, cereals, pancakes, etc.)"
[1] "Eggs"
[1] "Processed meat (other red or white meat such as lunch meat, ham, salami, bologna"
[1] "Red meat (beef, hamburger, pork, lamb)"
[1] "White meat (chicken, turkey, etc.)"
[1] "Shellfish (shrimp, lobster, scallops, etc.)"
[1] "Fish (fish nuggets, breaded fish, fish cakes, salmon, tuna, etc.)"
[1] "Sweets (pies, jam, chocolate, cake, cookies, etc.)"
r <- rbindlist(results)
ta <- rbindlist(tables)
sig <- r[p_abundance < 0.05/length(results)]
ggplot(sig) +
  aes(x=group, y=relative) +
  geom_jitter(data=ta[item %in% sig$item], 
              aes(y=log10(relative*1e6)), stroke=0, size=0.5, width=0.3) +
  geom_pointrange(aes(ymin=relative-sd, ymax=relative+sd), shape=21, fill="white") +
  geom_text(data=unique(sig, by="item"), aes(x=0, y=Inf, label=sprintf("p[lm]=%.2g\np[logit]=%.2g", p_abundance, p_prevalence)), size=3, vjust=1, hjust=0, nudge_y=0.5) +
  facet_wrap(~ item, nrow=1, labeller = labeller(item = label_wrap_gen(40))) +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1)) +
  labs(x="", y=expression(paste("Hibiscus abundance [",log[10],RPM, "]"))) +
  guides(color=F)

ggsave("figures/hibiscus.pdf", width=8, height=6)

Figure Data

openxlsx::write.xlsx(
  list(`Fig. ED6 (individual)` = ta, `Fig. ED6 (summaries)` = r),
  "figure_data/figED6.xlsx"
)

fwrite(metadata, "figure_data/ihmp_metadata.csv")
LS0tCnRpdGxlOiAiaUhNUCIKb3V0cHV0OiBodG1sX25vdGVib29rCi0tLQoKYGBge3J9CmxpYnJhcnkoZGF0YS50YWJsZSkKbGlicmFyeShtYWdyaXR0cikKbGlicmFyeShnZ3Bsb3QyKQp0aGVtZV9zZXQodGhlbWVfbWluaW1hbCgpKQoKbWV0YWRhdGEgPC0gZnJlYWQoImlobXAvZGF0YS9pYmRtZGJfaGVhbHRoeS5jc3YiKQpzcmFfZGF0YSA8LSBmcmVhZCgiaWhtcC9kYXRhL2liZG1kX2hlYWx0aHlfcnVudGFibGUudHN2IilbLCAuKFJ1biwgYExpYnJhcnkgTmFtZWApXQpzcmFfZGF0YVssICJFeHRlcm5hbCBJRCIgOj0gdHN0cnNwbGl0KGBMaWJyYXJ5IE5hbWVgLCAiX01HWCIpW1sxXV1dCm1ldGFkYXRhIDwtIG1ldGFkYXRhW3NyYV9kYXRhLCBvbj0iRXh0ZXJuYWwgSUQiLCBub21hdGNoPTBdCm1ldGFkYXRhCmBgYAoKYGBge3J9Cm1hbiA8LSBmcmVhZCgiZGIvZGF0YS9tYW5pZmVzdC5jc3YiKVssIC4oZGIsIHJhbmssIG1hdGNoZWRfdGF4aWQsIHNlcWxlbmd0aCldCgpmb29kcyA8LSBmcmVhZCgiaWhtcC9kYXRhL2Zvb2RfYWJ1bmRhbmNlLmNzdiIpCmZvb2RzIDwtIG1hbltmb29kcywgb249Im1hdGNoZWRfdGF4aWQiXQpmb29kc1ssICJ0cG0iIDo9IDEuMGU2ICogcmVhZHMgLyBhcy5kb3VibGUoc2VxbGVuZ3RoKV0KZm9vZHNbaXMubmEodHBtKSwgdHBtIDo9IDAuMF0KCmNvbnRlbnRzIDwtIGZyZWFkKCJpaG1wL2RhdGEvZm9vZF9jb250ZW50LmNzdiIpCmNvbnRlbnRzCmBgYAoKYGBge3J9CnRlc3RzIDwtIGxpc3QoCiAgdmVnZXRhYmxlcyA9IGxpc3QoZ3JvdXA9IlZlZ2V0YWJsZXMgKHNhbGFkLCB0b21hdG9lcywgb25pb25zLCBncmVlbnMsIGNhcnJvdHMsIHBlcHBlcnMsIGdyZWVuIGJlYW5zLCBldGMpIiwgZXg9ZXhwcmVzc2lvbihmb29kX2dyb3VwID09ICJWZWdldGFibGVzIikpLAogIGZydWl0cyA9IGxpc3QoZ3JvdXA9IkZydWl0cyAobm8ganVpY2UpIChBcHBsZXMsIHJhaXNpbnMsIGJhbmFuYXMsIG9yYW5nZXMsIHN0cmF3YmVycmllcywgYmx1ZWJlcnJpZXMiLCBleD1leHByZXNzaW9uKGZvb2RfZ3JvdXAgPT0gIkZydWl0cyIpKSwKICBiZWFucyA9IGxpc3QoZ3JvdXA9IkJlYW5zICh0b2Z1LCBzb3ksIHNveSBidXJnZXJzLCBsZW50aWxzLCBNZXhpY2FuIGJlYW5zLCBsaW1hIGJlYW5zIGV0YykiLCBleD1leHByZXNzaW9uKGZvb2Rfc3ViZ3JvdXAgPT0gIkJlYW5zIikpLAogIGB3aGl0ZSBtZWF0YCA9IGxpc3QoZ3JvdXA9IldoaXRlIG1lYXQgKGNoaWNrZW4sIHR1cmtleSwgZXRjLikiLCBleD1leHByZXNzaW9uKGZvb2Rfc3ViZ3JvdXAgPT0gIlBvdWx0cnkiKSksCiAgc2hlbGxmaXNoID0gbGlzdChncm91cD0iU2hlbGxmaXNoIChzaHJpbXAsIGxvYnN0ZXIsIHNjYWxsb3BzLCBldGMuKSIsIGV4PWV4cHJlc3Npb24oZm9vZF9zdWJncm91cCAlY2hpbiUgYygiTW9sbHVza3MiLCAiQ3J1c3RhY2VhbnMiKSkpLAogIGZpc2ggPSBsaXN0KGdyb3VwPSJGaXNoIChmaXNoIG51Z2dldHMsIGJyZWFkZWQgZmlzaCwgZmlzaCBjYWtlcywgc2FsbW9uLCB0dW5hLCBldGMuKSIsIGV4PWV4cHJlc3Npb24oZm9vZF9zdWJncm91cCA9PSAiRmlzaGVzIikpLAogIGByZWQgbWVhdGAgPSBsaXN0KGdyb3VwPSJSZWQgbWVhdCAoYmVlZiwgaGFtYnVyZ2VyLCBwb3JrLCBsYW1iKSIsIGV4PWV4cHJlc3Npb24oZm9vZF9zdWJncm91cCAlY2hpbiUgYygiU3dpbmUiLCAiQm92aW5lcyIsICJDYXByYWUiKSkpCikKYGBgCgoKCmBgYHtyLCBmaWcud2lkdGg9NCwgZmlnLmhlaWdodD0zfQpsaWJyYXJ5KHBhdGNod29yaykKCnRhYmxlcyA8LSBsaXN0KCkKCmZyZXFzIDwtIGMoCiAgYE5vLCBJIGRpZCBub3QgY29uc3VtZSB0aGVzZSBwcm9kdWN0cyBpbiB0aGUgbGFzdCA3IGRheXNgID0gMCwKICBgV2l0aGluIHRoZSBwYXN0IDQgdG8gNyBkYXlzYCA9IDEvNS41LAogIGBXaXRoaW4gdGhlIHBhc3QgMiB0byAzIGRheXNgID0gMS8yLjUsCiAgYFllc3RlcmRheSwgMSB0byAyIHRpbWVzYCA9IDEuNSwgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgIAogIGBZZXN0ZXJkYXksIDMgb3IgbW9yZSB0aW1lc2AgPSAzLAogIGBOQWAgPSAwCikKCnJlc3VsdHMgPC0gbGlzdCgpCnBsb3RzIDwtIGxpc3QoKQptZXJnZWQgPC0gbWV0YWRhdGFbZm9vZHMsIG9uPWMoYFJ1bmAgPSAic2FtcGxlX2lkIildCmZvciAoaSBpbiAxOmxlbmd0aCh0ZXN0cykpIHsKICBuYW1lIDwtIG5hbWVzKHRlc3RzKVtpXQogIHZhbHMgPC0gdGVzdHNbW2ldXQogIGR0IDwtIG1lcmdlZFssIC4oYWJ1bmRhbmNlPXN1bShyZWFkc1tldmFsKHZhbHMkZXgpXSksIGdyb3VwPS5TRFtbdmFscyRncm91cF1dWzFdLCB0b3RhbF9yZWFkcz10b3RhbF9yZWFkc1sxXSksIGJ5PSJFeHRlcm5hbCBJRCJdCiAgbnMgPC0gZHRbLCAuTiwgYnk9Imdyb3VwIl0gfD4gc2V0a2V5KGdyb3VwKQogIGR0IDwtIGR0W2dyb3VwICVjaGluJSBuYW1lcyhmcmVxcykgJiBuc1tncm91cCwgTl0gPj0gMTBdCiAgZHRbLCAiZnJlcXVlbmN5IiA6PSBmcmVxc1tncm91cF1dCiAgZHRbLCBncm91cCA6PSBmYWN0b3IoZ3JvdXAsIGxldmVscz1uYW1lcyhmcmVxcylbbmFtZXMoZnJlcXMpICVpbiUgZ3JvdXBdKV0KICBkdFssICJyZWxhdGl2ZSIgOj0gKGFidW5kYW5jZSArIDEpIC8gKHRvdGFsX3JlYWRzKV0KICBkdFssICJpdGVtIiA6PSBuYW1lXQogIHBsb3RzW1tuYW1lXV0gPC0gZ2dwbG90KGR0KSArIGFlcyh4PWdyb3VwLCB5PWxvZzEwKHJlbGF0aXZlKSkgKyAKICAgIGdlb21faml0dGVyKHdpZHRoPTAuMykgKyAKICAgIGdlb21fYm94cGxvdCh3aWR0aD0wLjEpICsgCiAgICBzdGF0X3Ntb290aChhZXMoZ3JvdXA9MSksIG1ldGhvZD0ibG0iKSArIAogICAgbGFicyh0aXRsZT1uYW1lLCB4PSJjb25zdW1wdGlvbiBmcmVxdWVuY3kgW3NlcnZpbmdzL2RheV0iLCB5PSJhYnVuZGFuY2UiKQogIHByaW50KG5hbWUpCiAgYWJtb2QgPC0gbG0obG9nMTAocmVsYXRpdmUpIH4gZnJlcXVlbmN5LCBkYXRhPWR0KQogIHByZW1vZCA8LSBnbG0oKGFidW5kYW5jZSA+IDApIH4gZnJlcXVlbmN5LCBkYXRhPWR0LCBmYW1pbHk9Ymlub21pYWwobGluaz0ibG9naXQiKSkKICB0YWJsZXNbW25hbWVdXSA8LSBkdAogIHJlc3VsdHNbW25hbWVdXSA8LSBkdFssIC4obj0uTiwgcmVsYXRpdmU9bWVhbihsb2cxMChyZWxhdGl2ZSoxZTYpKSwgc2Q9c2QobG9nMTAocmVsYXRpdmUqMWU2KSksIGRldGVjdGVkPXN1bShhYnVuZGFuY2UgPiAwKSwgaXRlbT1uYW1lKSwgYnk9Imdyb3VwIl0KICByZXN1bHRzW1tuYW1lXV1bLCAicF9hYnVuZGFuY2UiIDo9IHN1bW1hcnkoYWJtb2QpJGNvZWZmaWNpZW50c1syLCA0XV0KICByZXN1bHRzW1tuYW1lXV1bLCAicF9wcmV2YWxlbmNlIiA6PSBzdW1tYXJ5KHByZW1vZCkkY29lZmZpY2llbnRzWzIsIDRdXQp9CmBgYAoKYGBge3IsIGZpZy53aWR0aD0yLjUsIGZpZy5oZWlnaHQ9OX0KciA8LSByYmluZGxpc3QocmVzdWx0cykKdGEgPC0gcmJpbmRsaXN0KHRhYmxlcykKZ2dwbG90KHIpICsKICBhZXMoeD1ncm91cCwgeT1yZWxhdGl2ZSkgKwogIGdlb21faml0dGVyKGRhdGE9dGEsIGFlcyh5PWxvZzEwKHJlbGF0aXZlKjFlNikpLCBzdHJva2U9MCwgc2l6ZT0wLjUsIHdpZHRoPTAuMykgKwogIGdlb21fcG9pbnRyYW5nZShhZXMoeW1pbj1yZWxhdGl2ZS1zZCwgeW1heD1yZWxhdGl2ZStzZCksIHNoYXBlPTIxLCBmaWxsPSJ3aGl0ZSIpICsKICBnZW9tX3RleHQoZGF0YT11bmlxdWUociwgYnk9Iml0ZW0iKSwgYWVzKHg9MCwgeT1JbmYsIGxhYmVsPXNwcmludGYoInBbbG1dPSUuMmdcbnBbbG9naXRdPSUuMmciLCBwX2FidW5kYW5jZSwgcF9wcmV2YWxlbmNlKSksIHNpemU9Mywgdmp1c3Q9MSwgaGp1c3Q9MCwgbnVkZ2VfeT0wLjUpICsKICBmYWNldF93cmFwKH4gaXRlbSwgbmNvbD0yKSArCiAgdGhlbWUoYXhpcy50ZXh0LnggPSBlbGVtZW50X3RleHQoYW5nbGUgPSA5MCwgdmp1c3QgPSAwLjUsIGhqdXN0PTEpKSArCiAgbGFicyh4PSIiLCB5PWV4cHJlc3Npb24ocGFzdGUoImxvZyByZWxhdGl2ZSBhYnVuZGFuY2UgWyIsbG9nWzEwXSxSUE0sICJdIikpKSArCiAgZ3VpZGVzKGNvbG9yPUYpCmdnc2F2ZSgiZmlndXJlcy9paG1wX2Zvb2RzLnBkZiIsIHdpZHRoPTIuNSwgaGVpZ2h0PTkpCmBgYAoKCmBgYHtyLCBmaWcud2lkdGg9NSwgZmlnLmhlaWdodD0zfQpnZ3Bsb3QoZm9vZHNbc2VxbGVuZ3RoPjFlNl0sIGFlcyh4PXNlcWxlbmd0aC8xZTksIHk9cmVsYXRpdmVfYWJ1bmRhbmNlKSkgKyBnZW9tX3BvaW50KCkgKyBzdGF0X3Ntb290aChtZXRob2Q9ImxtIiwgY29sb3I9InRvbWF0byIpICsKICBsYWJzKHg9ImhhcGxvaWQgZ2Vub21lIHNpemUgW0dicHNdIiwgeT0icmVsYXRpdmUgZm9vZCBhYnVuZGFuY2UiKSArIHNjYWxlX3hfbG9nMTAoKSArIHNjYWxlX3lfbG9nMTAoKQpmb29kc1tzZXFsZW5ndGg+MWU2LCBjb3IudGVzdChyZWFkcywgc2VxbGVuZ3RoKV0KZ2dzYXZlKCJmaWd1cmVzL2dlbm9tZV9zaXplLnBkZiIsIHdpZHRoPTUsIGhlaWdodD0zKQpgYGAgCgojIyBEaWV0IGRpc3RhbmNlcwoKYGBge3J9CnNwZWNpZXMgPC0gZnJlYWQoImlobXAvZGF0YS9TX2NvdW50cy5jc3YiKVtncmVwbCgiQmFjdGVyaWF8QXJjaGFlYSIsIGxpbmVhZ2UpXQptaWNyb2Jpb21lX21hdCA8LSBkY2FzdChzcGVjaWVzLCBzYW1wbGVfaWQgfiBsaW5lYWdlLCB2YWx1ZS52YXIgPSAibmV3X2VzdF9yZWFkcyIsIGZpbGw9MCwgZnVuLmFnZ3JlZ2F0ZSA9IHN1bSkKc2lkcyA8LSBtaWNyb2Jpb21lX21hdFssIHRzdHJzcGxpdChzYW1wbGVfaWQsICJTXyIpW1syXV1dCm1pY3JvYmlvbWVfbWF0IDwtIGFzLm1hdHJpeChtaWNyb2Jpb21lX21hdFssICJzYW1wbGVfaWQiIDo9IE5VTExdKQpyb3duYW1lcyhtaWNyb2Jpb21lX21hdCkgPC0gc2lkcwpnb29kIDwtIGNvbE1lYW5zKG1pY3JvYmlvbWVfbWF0KSA+IDEwCm1pY3JvYmlvbWVfbWF0IDwtIG1pY3JvYmlvbWVfbWF0WywgZ29vZF0KIyBtaWNyb19yYXJlIDwtIHBoeWxvc2VxKG90dV90YWJsZShtaWNyb2Jpb21lX21hdCwgdGF4YV9hcmVfcm93cyA9IEYpKSB8PiByYXJlZnlfZXZlbl9kZXB0aCgxMDAwMDApIHw+IG90dV90YWJsZSgpCm1pY3JvYmlvbWVfcmVsYXRpdmUgPC0gbWljcm9iaW9tZV9tYXQgLyByb3dTdW1zKG1pY3JvYmlvbWVfbWF0KQpgYGAKCmBgYHtyfQpkaWV0IDwtIG1ldGFkYXRhWywgbmFtZXMobWV0YWRhdGEpWzcyOjkyXVstMTFdLCB3aXRoPUZdCmRpZXQgPC0gYXBwbHkoZGlldCwgMiwgZnVuY3Rpb24oeCkgaWZlbHNlKHggPT0gIiIsIDAsIGZyZXFzW3hdKSkKcm93bmFtZXMoZGlldCkgPC0gbWV0YWRhdGFbWyJSdW4iXV0KZGlldCA8LSBkaWV0W3Jvd1N1bXMoZGlldCkgPiAwLCBdCmRpZXQgPC0gZGlldFssIGNvbE1lYW5zKGRpZXQgPjApID4gMC4xXQpkaW0oZGlldCkKYGBgCgpgYGB7cn0KZm9vZGFiIDwtIGZvb2RzCmZvb2RhYlssICJpZCIgOj0gcGFzdGUobWF0Y2hlZF90YXhpZCwgc3BlY2llcywgd2lraXBlZGlhX2lkLCBzZXA9InwiKV0KZm9vZF9tYXQgPC0gZGNhc3QoZm9vZGFiLCBzYW1wbGVfaWQgfiBpZCwgdmFsdWUudmFyID0gInJlYWRzIiwgZmlsbD0wLCBmdW4uYWdncmVnYXRlID0gc3VtKQpzaWRzIDwtIGZvb2RfbWF0Wywgc2FtcGxlX2lkXQpmb29kX21hdCA8LSBhcy5tYXRyaXgoZm9vZF9tYXRbLCBzYW1wbGVfaWQgOj0gTlVMTF0pCnJvd25hbWVzKGZvb2RfbWF0KSA8LSBzaWRzCmZvb2RfcmVsYXRpdmUgPC0gZm9vZF9tYXRbcm93U3Vtcyhmb29kX21hdCkgPiAwLCBdCmZvb2RfcmVsYXRpdmUgPC0gZm9vZF9yZWxhdGl2ZSAvIHJvd1N1bXMoZm9vZF9yZWxhdGl2ZSkKZm9vZF9yZWxhdGl2ZSA8LSBmb29kX3JlbGF0aXZlW3Jvd1N1bXMoZm9vZF9yZWxhdGl2ZSkgPiAwLCBdCmZvb2RfcmVsYXRpdmUgPC0gZm9vZF9yZWxhdGl2ZVssIGNvbE1lYW5zKGZvb2RfcmVsYXRpdmUgPjApID4gMC4xXQpkaW0oZm9vZF9yZWxhdGl2ZSkKYGBgCgojIyBCZXRhIGRpdmVyc2l0eSB0ZXN0cwoKTGV0J3Mgd3JpdGUgYSBsaXR0bGUgZnVuY3Rpb24gdGhhdCBydW5zIHRoZSB0ZXN0IGZvciBtaWNyb2Jpb21lIDwtPiBvdGhlciBjb21wYXJpc29uIGFuZApwbG90cyBhbmQgcmV0dXJucyB0aGUgcmVzdWx0cy4KCmBgYHtyfQpsaWJyYXJ5KHZlZ2FuKQoKbWljcm9fbWFudGVsIDwtIGZ1bmN0aW9uKGZpcnN0LCBvdGhlciwgZmlyc3RuYW1lLCBvdGhlcm5hbWUpIHsKICBzaWRzIDwtIGludGVyc2VjdChyb3duYW1lcyhmaXJzdCksIHJvd25hbWVzKG90aGVyKSkKICBmaXJzdF9kaXN0IDwtIHZlZ2Rpc3QoZmlyc3Rbc2lkcywgXSwgImJyYXkiKQogIG90aGVyX2Rpc3QgPC0gdmVnZGlzdChvdGhlcltzaWRzLCBdLCAiYnJheSIpCiAgCiAgdGVzdCA8LSBtYW50ZWwoZmlyc3RfZGlzdCwgb3RoZXJfZGlzdCwgbWV0aG9kPSJwZWFyc29uIikKICByZXN1bHRzIDwtIGRhdGEudGFibGUoZmlyc3Q9YXMubnVtZXJpYyhmaXJzdF9kaXN0KSwgb3RoZXI9YXMubnVtZXJpYyhvdGhlcl9kaXN0KSkKICBwbCA8LSBnZ3Bsb3QocmVzdWx0cykgKwogICAgYWVzKHg9Zmlyc3QsIHk9b3RoZXIpICsKICAgIGdlb21fcG9pbnQoc2l6ZT0xLCBhbHBoYT0wLjI1LCBzdHJva2U9MCkgKwogICAgc3RhdF9zbW9vdGgobWV0aG9kPSJsbSIpICsKICAgIGxhYnMoeD1wYXN0ZShmaXJzdG5hbWUsICJkaXN0YW5jZSBbQnJheV0iKSwgeT1wYXN0ZShvdGhlcm5hbWUsICJkaXN0YW5jZSBbQnJheV0iKSkKICBwcmludChwbCkKICByZXR1cm4oZGF0YS50YWJsZSh4PWZpcnN0bmFtZSwgeT1vdGhlcm5hbWUsIHI9dGVzdCRzdGF0aXN0aWMsIHA9dGVzdCRzaWduaWYsIHBlcm1zPXRlc3RzJHBlcm11dGF0aW9ucykpCn0KYGBgCgpOb3cgd2UgcnVuIGl0IGZvciB0aGUgbWVhc3VyZXMuCgpgYGB7cn0KbWFudGVsX3Rlc3RzIDwtIGxpc3QoCiAgbWljcm9fbWFudGVsKG1pY3JvYmlvbWVfcmVsYXRpdmUsIGRpZXQsICJtaWNyb2Jpb21lIiwgIkZGUSIpLAogIG1pY3JvX21hbnRlbChtaWNyb2Jpb21lX3JlbGF0aXZlLCBmb29kX3JlbGF0aXZlLCAibWljcm9iaW9tZSIsICJNRURJIiksCiAgbWljcm9fbWFudGVsKGRpZXQsIGZvb2RfcmVsYXRpdmUsICJGRlEiLCAiTUVESSIpCikgfD4gcmJpbmRsaXN0KCkKbWFudGVsX3Rlc3RzCmBgYAoKIyMgRmlndXJlIGRhdGEKCmBgYHtyfQpmd3JpdGUodGEsICJmaWd1cmVfZGF0YS9pSE1QX3BvaW50cy5jc3YiKQpmd3JpdGUociwgImZpZ3VyZV9kYXRhL2lITVBfc3VtbWFyaWVzLmNzdiIpCmZ3cml0ZShtYW50ZWxfdGVzdHMsICJmaWd1cmVfZGF0YS9pSE1QX21hbnRlbC5jc3YiKQoKZndyaXRlKGZvb2RzW3NlcWxlbmd0aD4xZTZdLCAiZmlndXJlX2RhdGEvaWhtcF9nZW5vbWVfc2l6ZS5jc3YiKQpgYGAKCiMjIEFuYWx5c2lzIGZvciBoaWJpc2N1cyByZWFkcwoKYGBge3IsIGZpZy53aWR0aD00LCBmaWcuaGVpZ2h0PTN9CmxpYnJhcnkocGF0Y2h3b3JrKQoKdGFibGVzIDwtIGxpc3QoKQoKZnJlcXMgPC0gYygKICBgTm8sIEkgZGlkIG5vdCBjb25zdW1lIHRoZXNlIHByb2R1Y3RzIGluIHRoZSBsYXN0IDcgZGF5c2AgPSAwLAogIGBXaXRoaW4gdGhlIHBhc3QgNCB0byA3IGRheXNgID0gMS81LjUsCiAgYFdpdGhpbiB0aGUgcGFzdCAyIHRvIDMgZGF5c2AgPSAxLzIuNSwKICBgWWVzdGVyZGF5LCAxIHRvIDIgdGltZXNgID0gMS41LCAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgCiAgYFllc3RlcmRheSwgMyBvciBtb3JlIHRpbWVzYCA9IDMsCiAgYE5BYCA9IDAKKQp0ZXN0cyA8LSBuYW1lcyhtZXRhZGF0YSlbNzI6OTJdWy0xMV0KCnJlc3VsdHMgPC0gbGlzdCgpCnBsb3RzIDwtIGxpc3QoKQptZXJnZWQgPC0gbWV0YWRhdGFbZm9vZHMsIG9uPWMoYFJ1bmAgPSAic2FtcGxlX2lkIildCmZvciAoaSBpbiAxOmxlbmd0aCh0ZXN0cykpIHsKICBuYW1lIDwtIHRlc3RzW2ldCiAgZHQgPC0gbWVyZ2VkWywgLihhYnVuZGFuY2U9c3VtKHJlYWRzW2dlbnVzPT0iSGliaXNjdXMiXSksIGdyb3VwPS5TRFtbbmFtZV1dWzFdLCB0b3RhbF9yZWFkcz10b3RhbF9yZWFkc1sxXSksIGJ5PSJFeHRlcm5hbCBJRCJdCiAgbnMgPC0gZHRbLCAuTiwgYnk9Imdyb3VwIl0gfD4gc2V0a2V5KGdyb3VwKQogIGR0IDwtIGR0W2dyb3VwICVjaGluJSBuYW1lcyhmcmVxcykgJiBuc1tncm91cCwgTl0gPj0gMTBdCiAgZHRbLCAiZnJlcXVlbmN5IiA6PSBmcmVxc1tncm91cF1dCiAgZHRbLCBncm91cCA6PSBmYWN0b3IoZ3JvdXAsIGxldmVscz1uYW1lcyhmcmVxcylbbmFtZXMoZnJlcXMpICVpbiUgZ3JvdXBdKV0KICBkdFssICJyZWxhdGl2ZSIgOj0gKGFidW5kYW5jZSArIDEpIC8gKHRvdGFsX3JlYWRzKV0KICBkdFssICJpdGVtIiA6PSBuYW1lXQogIGlmIChkdFssIHVuaXF1ZU4oZ3JvdXApID09IDFdKSAgewogICAgbmV4dAogIH0KICBwbG90c1tbbmFtZV1dIDwtIGdncGxvdChkdCkgKyBhZXMoeD1ncm91cCwgeT1sb2cxMChyZWxhdGl2ZSkpICsgCiAgICBnZW9tX2ppdHRlcih3aWR0aD0wLjMpICsgCiAgICBnZW9tX2JveHBsb3Qod2lkdGg9MC4xKSArIAogICAgc3RhdF9zbW9vdGgoYWVzKGdyb3VwPTEpLCBtZXRob2Q9ImxtIikgKyAKICAgIGxhYnModGl0bGU9bmFtZSwgeD0iY29uc3VtcHRpb24gZnJlcXVlbmN5IFtzZXJ2aW5ncy9kYXldIiwgeT0iYWJ1bmRhbmNlIikKICBwcmludChuYW1lKQogIGFibW9kIDwtIGxtKGxvZzEwKHJlbGF0aXZlKSB+IGZyZXF1ZW5jeSwgZGF0YT1kdCkKICBwcmVtb2QgPC0gZ2xtKChhYnVuZGFuY2UgPiAwKSB+IGZyZXF1ZW5jeSwgZGF0YT1kdCwgZmFtaWx5PWJpbm9taWFsKGxpbms9ImxvZ2l0IikpCiAgdGFibGVzW1tuYW1lXV0gPC0gZHQKICByZXN1bHRzW1tuYW1lXV0gPC0gZHRbLCAuKG49Lk4sIHJlbGF0aXZlPW1lYW4obG9nMTAocmVsYXRpdmUqMWU2KSksIHNkPXNkKGxvZzEwKHJlbGF0aXZlKjFlNikpLCBkZXRlY3RlZD1zdW0oYWJ1bmRhbmNlID4gMCksIGl0ZW09bmFtZSksIGJ5PSJncm91cCJdCiAgcmVzdWx0c1tbbmFtZV1dWywgInBfYWJ1bmRhbmNlIiA6PSBzdW1tYXJ5KGFibW9kKSRjb2VmZmljaWVudHNbMiwgNF1dCiAgcmVzdWx0c1tbbmFtZV1dWywgInBfcHJldmFsZW5jZSIgOj0gc3VtbWFyeShwcmVtb2QpJGNvZWZmaWNpZW50c1syLCA0XV0KfQpgYGAKCmBgYHtyLCBmaWcud2lkdGg9OCwgZmlnLmhlaWdodD02fQpyIDwtIHJiaW5kbGlzdChyZXN1bHRzKQp0YSA8LSByYmluZGxpc3QodGFibGVzKQpzaWcgPC0gcltwX2FidW5kYW5jZSA8IDAuMDUvbGVuZ3RoKHJlc3VsdHMpXQpnZ3Bsb3Qoc2lnKSArCiAgYWVzKHg9Z3JvdXAsIHk9cmVsYXRpdmUpICsKICBnZW9tX2ppdHRlcihkYXRhPXRhW2l0ZW0gJWluJSBzaWckaXRlbV0sIAogICAgICAgICAgICAgIGFlcyh5PWxvZzEwKHJlbGF0aXZlKjFlNikpLCBzdHJva2U9MCwgc2l6ZT0wLjUsIHdpZHRoPTAuMykgKwogIGdlb21fcG9pbnRyYW5nZShhZXMoeW1pbj1yZWxhdGl2ZS1zZCwgeW1heD1yZWxhdGl2ZStzZCksIHNoYXBlPTIxLCBmaWxsPSJ3aGl0ZSIpICsKICBnZW9tX3RleHQoZGF0YT11bmlxdWUoc2lnLCBieT0iaXRlbSIpLCBhZXMoeD0wLCB5PUluZiwgbGFiZWw9c3ByaW50ZigicFtsbV09JS4yZ1xucFtsb2dpdF09JS4yZyIsIHBfYWJ1bmRhbmNlLCBwX3ByZXZhbGVuY2UpKSwgc2l6ZT0zLCB2anVzdD0xLCBoanVzdD0wLCBudWRnZV95PTAuNSkgKwogIGZhY2V0X3dyYXAofiBpdGVtLCBucm93PTEsIGxhYmVsbGVyID0gbGFiZWxsZXIoaXRlbSA9IGxhYmVsX3dyYXBfZ2VuKDQwKSkpICsKICB0aGVtZShheGlzLnRleHQueCA9IGVsZW1lbnRfdGV4dChhbmdsZSA9IDkwLCB2anVzdCA9IDAuNSwgaGp1c3Q9MSkpICsKICBsYWJzKHg9IiIsIHk9ZXhwcmVzc2lvbihwYXN0ZSgiSGliaXNjdXMgYWJ1bmRhbmNlIFsiLGxvZ1sxMF0sUlBNLCAiXSIpKSkgKwogIGd1aWRlcyhjb2xvcj1GKQpnZ3NhdmUoImZpZ3VyZXMvaGliaXNjdXMucGRmIiwgd2lkdGg9OCwgaGVpZ2h0PTYpCmBgYAoKIyMgRmlndXJlIERhdGEKCmBgYHtyfQpvcGVueGxzeDo6d3JpdGUueGxzeCgKICBsaXN0KGBGaWcuIEVENiAoaW5kaXZpZHVhbClgID0gdGEsIGBGaWcuIEVENiAoc3VtbWFyaWVzKWAgPSByKSwKICAiZmlndXJlX2RhdGEvZmlnRUQ2Lnhsc3giCikKCmZ3cml0ZShtZXRhZGF0YSwgImZpZ3VyZV9kYXRhL2lobXBfbWV0YWRhdGEuY3N2IikKYGBg